Load libraries

library(tigris)
library(zonebuilder)
library(tidyverse)
library(sf)
library(ggspatial)
library(units)

Read in data and create grids

boundary <- st_read("https://bostonopendata-boston.opendata.arcgis.com/datasets/142500a77e2a4dbeb94a86f7e0b568bc_9.geojson?outSR=%7B%22latestWkid%22%3A2249%2C%22wkid%22%3A102686%7D", quiet = TRUE) %>%
  st_transform(2249) # Transform to MA State plane

schools <- st_read("https://bostonopendata-boston.opendata.arcgis.com/datasets/1d9509a8b2fd485d9ad471ba2fdb1f90_0.geojson?outSR=%7B%22latestWkid%22%3A2249%2C%22wkid%22%3A102686%7D", quiet = TRUE) %>%
  st_transform(2249) # Transform to MA State plane

tracts <- tracts(state = "MA", county = "Suffolk") %>%
  st_transform(2249) # Transform to MA State plane

tracts <- tracts[boundary,]
  
grid <- st_sf(st_make_grid(boundary, n = c(20,20))) # Create a grid over Boston
grid <- grid[boundary,] # Filter for the cells that cover part of Boston

# Create a clockboard over Boston
clock <- zb_zone("Boston", distance = 0.5, distance_growth = 0, n_circles = 20) %>% 
  st_transform(2249) # Transform to MA State plane
clock <- clock[boundary,]  # Filter for  the cells that cover part of Boston

Show disagggregated points

plot1 <- ggplot(schools) +
  annotation_map_tile(zoomin = 0, progress = "none", type = "cartolight")  +
  geom_sf(size = 0.5, color = "darkred") +
  labs(caption = "Map tiles and data by OpenStreetMap") +
  theme_void()

plot1

jpeg("wk3_plot1.jpg", width = 7, height = 7, units = "in", res = 300)
plot1
dev.off()

Show Tract boundaries

plot2 <- ggplot(schools) +
  annotation_map_tile(zoomin = 0, progress = "none", type = "cartolight")  +
  geom_sf(size = 0.5, color = "darkred") +
  geom_sf(data = tracts, fill = NA) +
  labs(caption = "Map tiles and data by OpenStreetMap") +
  theme_void()

plot2

jpeg("wk3_plot2.jpg", width = 7, height = 7, units = "in", res = 300)
plot2
dev.off()

Aggregate to tracts

tracts <- tracts %>%
  mutate(num_schools = lengths(st_covers(tracts, schools))) %>%
  mutate(area = set_units(st_area(tracts), km^2)) %>%
  mutate(school_dens = as.numeric(num_schools / area))


plot3 <- ggplot(tracts) +
  annotation_map_tile(zoomin = 0, progress = "none", type = "cartolight")  +
  geom_sf(data = tracts, alpha = 0.5, aes(fill = school_dens)) +
  scale_fill_viridis_c(name = "Schools per square km") +
  labs(caption = "Map tiles and data by OpenStreetMap") +
  theme_void()

plot3

jpeg("wk3_plot3.jpg", width = 7, height = 7, units = "in", res = 300)
plot3
dev.off()

Show grid

plot4 <- ggplot(schools) +
  annotation_map_tile(zoomin = 0, progress = "none", type = "cartolight")  +
  geom_sf(size = 0.5, color = "darkred") +
  geom_sf(data = grid, fill = NA) +
  labs(caption = "Map tiles and data by OpenStreetMap") +
  theme_void()

plot4

jpeg("wk3_plot4.jpg", width = 7, height = 7, units = "in", res = 300)
plot4
dev.off()

Aggregate to grid

grid <- grid %>%
  mutate(num_schools = lengths(st_covers(grid, schools))) %>%
  mutate(area = set_units(st_area(grid), km^2)) %>%
  mutate(school_dens = as.numeric(num_schools / area))


plot5 <- ggplot(grid) +
  annotation_map_tile(zoomin = 0, progress = "none", type = "cartolight")  +
  geom_sf(data = grid, alpha = 0.5, aes(fill = school_dens)) +
  scale_fill_viridis_c(name = "Schools per square km") +
  labs(caption = "Map tiles and data by OpenStreetMap") +
  theme_void()

plot5

jpeg("wk3_plot5.jpg", width = 7, height = 7, units = "in", res = 300)
plot5
dev.off()

Show clockboard zones

plot6 <- ggplot(schools) +
  annotation_map_tile(zoomin = 0, progress = "none", type = "cartolight")  +
  geom_sf(size = 0.5, color = "darkred") +
  geom_sf(data = clock, fill = NA) +
  labs(caption = "Map tiles and data by OpenStreetMap") +
  theme_void()

plot6

jpeg("wk3_plot6.jpg", width = 7, height = 7, units = "in", res = 300)
plot6
dev.off()
## png 
##   2

Aggregate to clockboard zones

clock <- clock %>%
  mutate(num_schools = lengths(st_covers(clock, schools))) %>%
  mutate(area = set_units(st_area(clock), km^2)) %>%
  mutate(school_dens = as.numeric(num_schools / area))


plot7 <- ggplot(clock) +
  annotation_map_tile(zoomin = 0, progress = "none", type = "cartolight")  +
  geom_sf(alpha = 0.5, aes(fill = school_dens)) +
  scale_fill_viridis_c(name = "Schools per square km") +
  labs(caption = "Map tiles and data by OpenStreetMap") +
  theme_void()

plot7

jpeg("wk3_plot7.jpg", width = 7, height = 7, units = "in", res = 300)
plot7
dev.off()